function meanX = meanProj(theta,setup)

nx             = setup.nx;
ny             = setup.ny;
scaleX_in_g    = setup.scaleX_in_g;
g0             = theta(:,1);    %contains the steady state AND the risk-adj for the controls
% First order terms
gx             = theta(:,1+1:1+nx)./repmat(setup.scaleX_in_g',ny,1);

% Nonlinear terms
gxx           = zeros(ny,nx,nx);
idx           = 0;
if setup.appOrder > 1
    for i=1:nx
        for j=i:nx
            idx = idx + 1;
            gxx(:,i,j) = theta(:,1+nx+idx)/(scaleX_in_g(i,1)*scaleX_in_g(j,1));
            gxx(:,j,i) = gxx(:,i,j);
        end
    end
end

% The states
h0   = zeros(nx,1);
hx   = zeros(nx,nx);
hxx  = zeros(nx,nx,nx);
mx   = setup.mx;
myx  = setup.myx;
if setup.mx > 0
    % Pure Endogenous states
    error('projectStep.m: cannot handled mx > 0')
end
% recall states are in dev from steady state
if setup.myx > 0
    % Laggend controls as states
    h0(mx+1:mx+myx,1)       = g0(1:myx,1)- setup.gSS(1:myx,1);  %contains risk adjustment for the states 
    hx(mx+1:mx+myx,:)       = gx(1:myx,:);
    hxx(mx+1:mx+myx,:,:)    = gxx(1:myx,:,:);
end
% Shocks are linear
hx(setup.mx+setup.myx+1:nx,:) = setup.hx(mx+myx+1:nx,:); 

varX     = reshape((eye(nx^2)-kron(hx,hx))\reshape((setup.eta*setup.eta'),nx^2,1),nx,nx);
meanX    = (eye(nx)-hx)\(h0+reshape(hxx,nx,nx^2)*varX(:));


end

